WL2 Growth + Climate

Add climate to best base model

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(brms)
## Loading required package: Rcpp
## Loading 'brms' package (version 2.21.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
## 
## Attaching package: 'brms'
## 
## The following object is masked from 'package:stats':
## 
##     ar
library(zoo) #rollmean
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

Load the growth data

all_wl2_size_all_nobuffers <- read_csv("../output/WL2_Traits/WL2-2023_Size_Combined.csv")
## Rows: 17336 Columns: 10
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (5): block, Genotype, pop.mf, parent.pop, survey.notes
## dbl  (4): mf, rep, height.cm, long.leaf.cm
## date (1): survey_date
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Pop Location Info

wl2_gowers_2023 <- read_csv("../output/Climate/Gowers_WL2.csv") %>% 
  select(parent.pop:GrwSsn_GD, Wtr_Year_GD) %>% 
  pivot_wider(names_from = TimePd, values_from = c(GrwSsn_GD, Wtr_Year_GD)) 
## Rows: 46 Columns: 12
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): parent.pop, elevation.group, TimePd
## dbl (9): elev_m, Lat, Long, GrwSsn_GD, GrwSsn_FLINT_GD, GrwSsn_BIOCLIM_GD, W...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(wl2_gowers_2023)
## # A tibble: 6 × 9
##   parent.pop elevation.group elev_m   Lat  Long GrwSsn_GD_Recent
##   <chr>      <chr>            <dbl> <dbl> <dbl>            <dbl>
## 1 YO4        High             2158.  37.8 -120.            0.216
## 2 WL2        High             2020.  38.8 -120.            0.242
## 3 BH         Low               511.  37.4 -120.            0.264
## 4 SQ3        High             2373.  36.7 -119.            0.265
## 5 SQ2        Mid              1934.  36.7 -119.            0.270
## 6 SQ1        Mid              1921.  36.6 -119.            0.271
## # ℹ 3 more variables: GrwSsn_GD_Historical <dbl>, Wtr_Year_GD_Recent <dbl>,
## #   Wtr_Year_GD_Historical <dbl>

Merge

all_wl2_size_loc <- left_join(all_wl2_size_all_nobuffers, wl2_gowers_2023) %>% 
  mutate(day=if_else(survey_date=="2023-07-03" | survey_date=="2023-07-06", 0, #day 0 = pre-transplant
                 as.numeric(survey_date-min(survey_date)-7))) #days in the field starts on 7/11 (day all plants taken to the field) rather than the exact planting date 
## Joining with `by = join_by(parent.pop)`
head(all_wl2_size_loc)
## # A tibble: 6 × 19
##   survey_date block Genotype pop.mf parent.pop    mf   rep height.cm
##   <date>      <chr> <chr>    <chr>  <chr>      <dbl> <dbl>     <dbl>
## 1 2023-07-03  <NA>  CP2_1_1  CP2_1  CP2            1     1       0.5
## 2 2023-07-03  <NA>  CP2_1_2  CP2_1  CP2            1     2       0.7
## 3 2023-07-03  <NA>  CP2_1_3  CP2_1  CP2            1     3       1.1
## 4 2023-07-03  <NA>  CP2_1_4  CP2_1  CP2            1     4       0.8
## 5 2023-07-03  <NA>  CP2_1_5  CP2_1  CP2            1     5       0.9
## 6 2023-07-03  <NA>  CP2_1_6  CP2_1  CP2            1     6       1  
## # ℹ 11 more variables: long.leaf.cm <dbl>, survey.notes <chr>,
## #   elevation.group <chr>, elev_m <dbl>, Lat <dbl>, Long <dbl>,
## #   GrwSsn_GD_Recent <dbl>, GrwSsn_GD_Historical <dbl>,
## #   Wtr_Year_GD_Recent <dbl>, Wtr_Year_GD_Historical <dbl>, day <dbl>

Data set for the model (remove WV and shrinkage months)

all_wl2_size_for_models <- all_wl2_size_loc %>% 
  filter(parent.pop!="WV") %>% 
  group_by(Genotype) %>% 
  mutate(height_next = lead(height.cm, order_by = survey_date),
         height_diff=height_next-height.cm,
         shrink_day=if_else(height_diff>-2, NA, day),
         min_shrink_day=min(shrink_day, na.rm = TRUE)) %>% #deals with problem of if something shrinks and then starts to grow again
  filter(day<min_shrink_day | min_shrink_day==Inf) %>% 
  ungroup() %>% 
  drop_na(height.cm)
## Warning: There were 1392 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `min_shrink_day = min(shrink_day, na.rm = TRUE)`.
## ℹ In group 1: `Genotype = "BH_1_1"`.
## Caused by warning in `min()`:
## ! no non-missing arguments to min; returning Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1391 remaining warnings.

Read in soil temp

temp_WL2 <- read_csv("../input/WL2_Data/WL2_2022_2023_iButton_Data_Corrected.csv") %>%
  select(-`...3`) %>%
  mutate(Date_Time = mdy_hm(Date_Time)) %>%
  filter(Date_Time > ymd("2023-07-06"))
## New names:
## Rows: 14253 Columns: 4
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (3): Bed, Date_Time, ...3 dbl (1): SoilTemp
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...3`
head(temp_WL2)
## # A tibble: 6 × 3
##   Bed   Date_Time           SoilTemp
##   <chr> <dttm>                 <dbl>
## 1 A_1   2023-07-13 14:46:00     27.5
## 2 A_1   2023-07-13 15:46:00     32  
## 3 A_1   2023-07-13 16:46:00     31.5
## 4 A_1   2023-07-13 17:46:00     28.5
## 5 A_1   2023-07-13 18:46:00     25  
## 6 A_1   2023-07-13 19:46:00     22
skimr::skim(temp_WL2)
Data summary
Name temp_WL2
Number of rows 10240
Number of columns 3
_______________________
Column type frequency:
character 1
numeric 1
POSIXct 1
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
Bed 0 1 3 3 0 5 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
SoilTemp 0 1 16.22 8.51 1.5 10.5 14 20 50.5 ▅▇▂▁▁

Variable type: POSIXct

skim_variable n_missing complete_rate min max median n_unique
Date_Time 0 1 2023-07-13 14:40:00 2023-10-06 21:48:00 2023-08-25 06:14:00 10240

Average across beds.

temp_WL2_summary <- temp_WL2 %>%
  group_by(Date_Time) %>% 
  summarise(AvgSoilTemp=mean(SoilTemp, na.rm=TRUE)) %>% #avg across beds 
  mutate(Date=as.Date(Date_Time)) %>%
  filter(Date != min(Date), Date != max(Date)) %>% # trim ragged ends
  group_by(Date) %>% #summmarise hourly data by date
  summarize(
    min_temp_d = min(AvgSoilTemp),
    max_temp_d = max(AvgSoilTemp),
    mean_temp_d = mean(AvgSoilTemp)
  ) %>%
  mutate(
    across(ends_with("temp_d"), \(x) rollmean(x, k = 7, align = "right", fill = NA), .names="{.col}1_7"),
    across(ends_with("temp_d"), \(x) rollmean(x, k = 13, align = "right", fill = NA), .names="{.col}1_13") # 13 so I can get the first survey date
  ) 

temp_WL2_summary
## # A tibble: 84 × 10
##    Date       min_temp_d max_temp_d mean_temp_d min_temp_d1_7 max_temp_d1_7
##    <date>          <dbl>      <dbl>       <dbl>         <dbl>         <dbl>
##  1 2023-07-14        8.5       31          17.0          NA            NA  
##  2 2023-07-15       11         35          18.7          NA            NA  
##  3 2023-07-16       11         36          19.1          NA            NA  
##  4 2023-07-17       14         36          20.5          NA            NA  
##  5 2023-07-18       13         36          19.8          NA            NA  
##  6 2023-07-19       12.5       36          19.3          NA            NA  
##  7 2023-07-20       11.5       45          21.4          11.6          36.4
##  8 2023-07-21       11         47.5        21.9          12            38.8
##  9 2023-07-22       11.5       48.5        22.5          12.1          40.7
## 10 2023-07-23       14         48.5        23.4          12.5          42.5
## # ℹ 74 more rows
## # ℹ 4 more variables: mean_temp_d1_7 <dbl>, min_temp_d1_13 <dbl>,
## #   max_temp_d1_13 <dbl>, mean_temp_d1_13 <dbl>

Read in soil moisture

moisture_WL2 <- read_csv("../input/WL2_Data/WL2_2023_Bed_C_Soil_Moisture_Corrected.csv") %>%
  mutate(Date_Time = mdy_hm(Date_Time)) 
## Rows: 2209 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Date_Time
## dbl (5): Port_1, Port_2, Port_3, Port_4, Port_5
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(moisture_WL2)
## # A tibble: 6 × 6
##   Date_Time           Port_1 Port_2 Port_3 Port_4 Port_5
##   <dttm>               <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
## 1 2023-07-20 12:00:00  0.071  0.141  0.075  0.055  0.092
## 2 2023-07-20 13:00:00  0.07   0.142  0.077  0.056  0.094
## 3 2023-07-20 14:00:00  0.074  0.142  0.079  0.056  0.094
## 4 2023-07-20 15:00:00  0.076  0.141  0.081  0.055  0.094
## 5 2023-07-20 16:00:00  0.077  0.14   0.082  0.054  0.093
## 6 2023-07-20 17:00:00  0.077  0.136  0.081  0.052  0.091
skimr::skim(moisture_WL2)
Data summary
Name moisture_WL2
Number of rows 2209
Number of columns 6
_______________________
Column type frequency:
numeric 5
POSIXct 1
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Port_1 0 1 0.13 0.06 0.00 0.08 0.14 0.17 0.36 ▆▇▇▁▁
Port_2 0 1 0.10 0.02 0.04 0.08 0.10 0.12 0.16 ▁▅▇▇▂
Port_3 0 1 0.10 0.05 -0.01 0.06 0.12 0.13 0.21 ▃▂▅▇▁
Port_4 0 1 0.03 0.02 -0.05 0.01 0.03 0.05 0.10 ▁▂▇▆▁
Port_5 0 1 0.08 0.04 -0.03 0.06 0.09 0.10 0.19 ▂▂▇▅▁

Variable type: POSIXct

skim_variable n_missing complete_rate min max median n_unique
Date_Time 0 1 2023-07-20 12:00:00 2023-10-20 12:00:00 2023-09-04 12:00:00 2209

Will average across the ports

moisture_WL2_summary <- moisture_WL2 %>%
  mutate(Port_1=if_else(Port_1<0, 0, Port_1),
         Port_2=if_else(Port_2<0, 0, Port_2),
         Port_3=if_else(Port_3<0, 0, Port_3),
         Port_4=if_else(Port_4<0, 0, Port_4),
         Port_5=if_else(Port_5<0, 0, Port_5)) %>%  #convert negatives to zeros (negatives result from extremely dry soil --> air pockets)
  rowwise() %>%
  mutate(s_moisture = mean(c_across(-Date_Time)) ) %>%
  select(Date_Time, s_moisture) %>%
  mutate(Date=as.Date(Date_Time)) %>%
  group_by(Date) %>%
  summarize(
    mean_moisture_d = mean(s_moisture)
  ) %>%
  mutate(
    s_moisture_1_7 = rollmean(mean_moisture_d, k = 7, align = "right", fill = "extend"),
    s_moisture_1_13 = rollmean(mean_moisture_d, k = 13, align = "right", fill = "extend")
  )

moisture_WL2_summary
## # A tibble: 93 × 4
##    Date       mean_moisture_d s_moisture_1_7 s_moisture_1_13
##    <date>               <dbl>          <dbl>           <dbl>
##  1 2023-07-20          0.0853         0.0714          0.0571
##  2 2023-07-21          0.0809         0.0714          0.0571
##  3 2023-07-22          0.0754         0.0714          0.0571
##  4 2023-07-23          0.0711         0.0714          0.0571
##  5 2023-07-24          0.0668         0.0714          0.0571
##  6 2023-07-25          0.0622         0.0714          0.0571
##  7 2023-07-26          0.0581         0.0714          0.0571
##  8 2023-07-27          0.0537         0.0669          0.0571
##  9 2023-07-28          0.0484         0.0622          0.0571
## 10 2023-07-29          0.0436         0.0577          0.0571
## # ℹ 83 more rows

add climate data

all_wl2_size_for_models_climate <- all_wl2_size_for_models %>% 
  left_join(temp_WL2_summary, by = c("survey_date" = "Date")) %>%
  left_join(moisture_WL2_summary, by = c("survey_date" = "Date")) %>% 
  filter(day>0) #remove day pre-field size b/c no climate data 

all_wl2_size_for_models_climate %>% 
  ggplot(aes(group=Genotype, x=survey_date, y=height.cm, col=elev_m)) + 
  scale_colour_gradient(low = "#F5A540", high = "#0043F0") +
  geom_line() + facet_wrap(~parent.pop, scales = "free")

#without day 0, min is different...

Priors

prior1_wl2 <- c(set_prior("normal(5,3)", nlpar="Hmin"),
            set_prior("normal(60,15)",  nlpar="Hmax"),
            set_prior("gamma(2,2)", nlpar = "k", lb = 0), #gamma constrains to be > 0 
            set_prior("gamma(20,3)", nlpar = "delta", lb = 0)
)
            #first number = avg, second number = the space around that average to search 

Base Models

#best base model:
base_f1 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
                  Hmin + k + delta ~ 1,
                  Hmax ~ (1|parent.pop) + (1|Genotype),
                  nl=TRUE)

base_f2 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
                  k + delta ~ 1,
                  Hmin + Hmax ~ (1|parent.pop) + (1|Genotype),
                  nl=TRUE)

Base Model Fits

Base model 1 - Hmax vary by pop

base_fit1 <- brm(formula=base_f1, data=all_wl2_size_for_models_climate, prior=prior1_wl2, iter = 4000, cores=4,  
                 control = list(max_treedepth = 12))
## Compiling Stan program...
## Trying to compile a simple C file
## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG   -I"/usr/lib/R/site-library/Rcpp/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/unsupported"  -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/"  -I"/usr/lib/R/site-library/StanHeaders/include/"  -I"/usr/lib/R/site-library/RcppParallel/include/"  -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
##                  from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
##                  from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
##                  from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
##   679 | #include <cmath>
##       |          ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1
## Start sampling
saveRDS(base_fit1, file = "brm_object_wl2_base1_rptedmeas.rds")
summary(base_fit1)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta) 
##          Hmin ~ 1
##          k ~ 1
##          delta ~ 1
##          Hmax ~ (1 | parent.pop) + (1 | Genotype)
##    Data: all_wl2_size_for_models_climate (Number of observations: 6038) 
##   Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup draws = 8000
## 
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmax_Intercept)     7.33      1.49     5.39    11.22 1.00     3514     2939
## 
## ~parent.pop (Number of levels: 22) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmax_Intercept)    11.47      2.97     7.38    18.85 1.00     2458     2860
## 
## Regression Coefficients:
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Hmin_Intercept      2.64      0.07     2.50     2.77 1.00     8880     6398
## k_Intercept         0.41      0.17     0.13     0.77 1.00     3610     2851
## delta_Intercept     0.59      0.03     0.54     0.64 1.00     3702     3209
## Hmax_Intercept      9.32      2.93     4.49    15.99 1.00     1093     2238
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.07      0.01     1.05     1.09 1.00    10226     6189
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(base_fit1, nvariables = 3, ask=FALSE)

pred.df <- expand_grid(day=min(all_wl2_size_for_models_climate$day):max(all_wl2_size_for_models_climate$day),
                       Genotype=unique(all_wl2_size_for_models_climate$Genotype)) %>% 
  separate(Genotype, c("parent.pop", "mf","rep"), remove = FALSE)
pred.df2 <- all_wl2_size_for_models_climate %>% select(day, parent.pop, Genotype)
base_fit1.predictions <- cbind(pred.df2, prediction=predict(base_fit1, newdata = pred.df2)[,"Estimate"]) %>%
  full_join(all_wl2_size_for_models_climate, by=c("day", "parent.pop", "Genotype"))

base_fit1.predictions %>% ggplot(aes(x=day)) +
  geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
  geom_line(aes(group=Genotype, y=prediction, color=parent.pop)) +
  facet_wrap(~parent.pop, scales = "free")

ggsave("../output/WL2_GrowthClim_Predictions_Base1_Hmax_RptedMeas.png", width = 14, height = 8, units = "in")

Base model 2 - Hmin + Hmax vary by pop

base_fit2 <- brm(formula=base_f2, data=all_wl2_size_for_models_climate, prior=prior1_wl2, iter = 4000, cores=4,
                 control = list(max_treedepth = 12))
## Compiling Stan program...
## Trying to compile a simple C file
## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG   -I"/usr/lib/R/site-library/Rcpp/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/unsupported"  -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/"  -I"/usr/lib/R/site-library/StanHeaders/include/"  -I"/usr/lib/R/site-library/RcppParallel/include/"  -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
##                  from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
##                  from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
##                  from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
##   679 | #include <cmath>
##       |          ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1
## Start sampling
saveRDS(base_fit2, file = "brm_object_wl2_base2_rptedmeas.rds")
summary(base_fit2)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta) 
##          k ~ 1
##          delta ~ 1
##          Hmin ~ (1 | parent.pop) + (1 | Genotype)
##          Hmax ~ (1 | parent.pop) + (1 | Genotype)
##    Data: all_wl2_size_for_models_climate (Number of observations: 6038) 
##   Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup draws = 8000
## 
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept)     1.86      0.05     1.76     1.96 1.00     5146     6543
## sd(Hmax_Intercept)     5.19      0.41     4.49     6.08 1.00     2519     3973
## 
## ~parent.pop (Number of levels: 22) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept)     1.39      0.25     1.00     1.95 1.00     3936     5209
## sd(Hmax_Intercept)     8.40      1.43     6.12    11.67 1.00     4214     4876
## 
## Regression Coefficients:
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## k_Intercept         0.83      0.10     0.64     1.02 1.00     2489     4084
## delta_Intercept     1.03      0.05     0.95     1.13 1.00     2279     3882
## Hmin_Intercept      3.24      0.31     2.64     3.87 1.00     2850     4409
## Hmax_Intercept      7.01      1.85     3.48    10.83 1.00     1980     3037
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.84      0.01     0.82     0.86 1.00     7479     6494
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(base_fit2, nvariables = 3, ask=FALSE)

pred.df <- expand_grid(day=min(all_wl2_size_for_models_climate$day):max(all_wl2_size_for_models_climate$day),
                       Genotype=unique(all_wl2_size_for_models_climate$Genotype)) %>% 
  separate(Genotype, c("parent.pop", "mf","rep"), remove = FALSE)
pred.df2 <- all_wl2_size_for_models_climate %>% select(day, parent.pop, Genotype)
base_fit2.predictions <- cbind(pred.df2, prediction=predict(base_fit2, newdata = pred.df2)[,"Estimate"]) %>%
  full_join(all_wl2_size_for_models_climate, by=c("day", "parent.pop", "Genotype"))

base_fit2.predictions %>% ggplot(aes(x=day)) +
  geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
  geom_line(aes(group=Genotype, y=prediction, color=parent.pop)) +
  facet_wrap(~parent.pop, scales = "free")

ggsave("../output/WL2_GrowthClim_Predictions_Base2_HmaxHmin_RptedMeas.png", width = 14, height = 8, units = "in")

Model comparison

base_fit1 <- add_criterion(base_fit1, "loo")
## Warning: Found 93 observations with a pareto_k > 0.7 in model 'base_fit1'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
base_fit2 <- add_criterion(base_fit2, "loo")
## Warning: Found 308 observations with a pareto_k > 0.7 in model 'base_fit2'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
loo_compare(base_fit1, base_fit2)
##           elpd_diff se_diff
## base_fit2     0.0       0.0
## base_fit1 -1000.5     100.9
#model with Hmin is better 

Temp Models

temp_f2 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
                  delta ~ 1,
                  k ~ min_temp_d1_7,
                  Hmin + Hmax ~ (1|parent.pop)+ (1|Genotype),
                  nl=TRUE)

temp_f3 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
                  delta ~ 1,
                  k ~ mean_temp_d1_7,
                  Hmin + Hmax ~ (1|parent.pop) + (1|Genotype),
                  nl=TRUE)

temp_f4 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
                  delta ~ 1,
                  k ~ max_temp_d1_7,
                  Hmin + Hmax ~ (1|parent.pop) + (1|Genotype),
                  nl=TRUE)

temp_f5 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
                  delta ~ 1,
                  k ~ min_temp_d1_13, 
                  Hmin + Hmax ~ (1|parent.pop) + (1|Genotype),
                  nl=TRUE)

temp_f6 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
                  delta ~ 1,
                  k ~ mean_temp_d1_13,
                  Hmin + Hmax ~ (1|parent.pop) + (1|Genotype),
                  nl=TRUE)

temp_f7 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
                  delta ~ 1,
                  k ~ max_temp_d1_13,
                  Hmin + Hmax ~ (1|parent.pop) + (1|Genotype),
                  nl=TRUE)

Temp Model Fits

Min D1-7 (fit2)

all_wl2_size_for_models_temp_fit2 <- all_wl2_size_for_models_climate %>% drop_na(min_temp_d1_7)

temp_fit2 <- brm(formula=temp_f2, data=all_wl2_size_for_models_temp_fit2, prior=prior1_wl2, iter = 3800, cores=4)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG   -I"/usr/lib/R/site-library/Rcpp/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/unsupported"  -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/"  -I"/usr/lib/R/site-library/StanHeaders/include/"  -I"/usr/lib/R/site-library/RcppParallel/include/"  -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
##                  from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
##                  from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
##                  from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
##   679 | #include <cmath>
##       |          ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1
## Start sampling
#Warning: There were 257 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceededWarning: Examine the pairs() plot to diagnose sampling problems
summary(temp_fit2) #some ESS kind of low 
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta) 
##          delta ~ 1
##          k ~ min_temp_d1_7
##          Hmin ~ (1 | parent.pop) + (1 | Genotype)
##          Hmax ~ (1 | parent.pop) + (1 | Genotype)
##    Data: all_wl2_size_for_models_temp_fit2 (Number of observations: 5403) 
##   Draws: 4 chains, each with iter = 3800; warmup = 1900; thin = 1;
##          total post-warmup draws = 7600
## 
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept)     1.78      0.05     1.68     1.88 1.00     4557     6223
## sd(Hmax_Intercept)     6.98      0.95     5.50     9.14 1.00     2483     3873
## 
## ~parent.pop (Number of levels: 22) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept)     1.48      0.26     1.07     2.07 1.00     3818     5009
## sd(Hmax_Intercept)    10.60      2.12     7.28    15.54 1.00     3015     4361
## 
## Regression Coefficients:
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## delta_Intercept     1.08      0.05     0.99     1.17 1.00     2516     5128
## k_Intercept         0.62      0.10     0.43     0.82 1.00     2347     3975
## k_min_temp_d1_7     0.00      0.00     0.00     0.01 1.00     8888     5563
## Hmin_Intercept      3.28      0.32     2.65     3.91 1.00     2207     3628
## Hmax_Intercept      8.41      2.49     3.87    13.69 1.00     1339     2854
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.76      0.01     0.74     0.77 1.00     7138     6042
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(temp_fit2, nvariables = 3, ask=FALSE) #some Hmax chain weirdness 

Mean D1-7 (fit3)

all_wl2_size_for_models_temp_fit3 <- all_wl2_size_for_models_climate %>% drop_na(mean_temp_d1_7)

temp_fit3 <- brm(formula=temp_f3, data=all_wl2_size_for_models_temp_fit3, prior=prior1_wl2, iter = 3800, cores=4)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG   -I"/usr/lib/R/site-library/Rcpp/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/unsupported"  -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/"  -I"/usr/lib/R/site-library/StanHeaders/include/"  -I"/usr/lib/R/site-library/RcppParallel/include/"  -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
##                  from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
##                  from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
##                  from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
##   679 | #include <cmath>
##       |          ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1
## Start sampling
summary(temp_fit3)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta) 
##          delta ~ 1
##          k ~ mean_temp_d1_7
##          Hmin ~ (1 | parent.pop) + (1 | Genotype)
##          Hmax ~ (1 | parent.pop) + (1 | Genotype)
##    Data: all_wl2_size_for_models_temp_fit3 (Number of observations: 5403) 
##   Draws: 4 chains, each with iter = 3800; warmup = 1900; thin = 1;
##          total post-warmup draws = 7600
## 
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept)     1.77      0.05     1.68     1.87 1.00     4391     5105
## sd(Hmax_Intercept)     6.89      0.92     5.46     9.07 1.00     2181     3419
## 
## ~parent.pop (Number of levels: 22) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept)     1.47      0.25     1.06     2.05 1.00     3683     5031
## sd(Hmax_Intercept)    10.44      2.08     7.26    15.27 1.00     3310     4698
## 
## Regression Coefficients:
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## delta_Intercept      1.13      0.05     1.04     1.22 1.00     3105     5145
## k_Intercept          0.59      0.09     0.41     0.77 1.00     2143     3568
## k_mean_temp_d1_7     0.01      0.00     0.00     0.01 1.00     5490     4379
## Hmin_Intercept       3.26      0.32     2.64     3.89 1.00     2274     3285
## Hmax_Intercept       8.36      2.43     3.93    13.54 1.00     1521     2778
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.75      0.01     0.74     0.77 1.00     5932     6072
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(temp_fit3, nvariables = 3, ask=FALSE)

#temp_fit3.predictions <- cbind(pred.df, prediction=predict(temp_fit3, newdata = pred.df)[,"Estimate"]) %>%
#  full_join(all_wl2_size_for_models, by=c("day", "parent.pop"))
#
#temp_fit3.predictions %>% ggplot(aes(x=day)) +
#  geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
#  geom_line(aes(y=prediction, color=parent.pop)) +
#  facet_wrap(~parent.pop)

Max D1-7 (fit4)

all_wl2_size_for_models_temp_fit4 <- all_wl2_size_for_models_climate %>% drop_na(max_temp_d1_7)

temp_fit4 <- brm(formula=temp_f4, data=all_wl2_size_for_models_temp_fit4, prior=prior1_wl2, iter = 3800, cores=4)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG   -I"/usr/lib/R/site-library/Rcpp/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/unsupported"  -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/"  -I"/usr/lib/R/site-library/StanHeaders/include/"  -I"/usr/lib/R/site-library/RcppParallel/include/"  -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
##                  from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
##                  from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
##                  from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
##   679 | #include <cmath>
##       |          ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1
## Start sampling
summary(temp_fit4)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta) 
##          delta ~ 1
##          k ~ max_temp_d1_7
##          Hmin ~ (1 | parent.pop) + (1 | Genotype)
##          Hmax ~ (1 | parent.pop) + (1 | Genotype)
##    Data: all_wl2_size_for_models_temp_fit4 (Number of observations: 5403) 
##   Draws: 4 chains, each with iter = 3800; warmup = 1900; thin = 1;
##          total post-warmup draws = 7600
## 
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept)     1.80      0.05     1.71     1.90 1.00     4281     5523
## sd(Hmax_Intercept)     5.63      0.64     4.62     7.12 1.00     2106     3038
## 
## ~parent.pop (Number of levels: 22) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept)     1.37      0.24     0.99     1.91 1.00     4333     5168
## sd(Hmax_Intercept)     8.74      1.65     6.14    12.60 1.00     3431     4947
## 
## Regression Coefficients:
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## delta_Intercept     1.06      0.04     0.98     1.14 1.00     2880     4393
## k_Intercept         0.74      0.11     0.53     0.95 1.00     2008     3185
## k_max_temp_d1_7     0.00      0.00     0.00     0.00 1.00     2687     4353
## Hmin_Intercept      3.20      0.31     2.59     3.80 1.00     2968     4317
## Hmax_Intercept      7.32      1.99     3.67    11.46 1.00     1788     3177
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.75      0.01     0.74     0.77 1.00     5457     5091
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(temp_fit4, nvariables = 3, ask=FALSE)

#temp_fit4.predictions <- cbind(pred.df, prediction=predict(temp_fit4, newdata = pred.df)[,"Estimate"]) %>%
#  full_join(all_wl2_size_for_models, by=c("day", "parent.pop"))
#
#temp_fit4.predictions %>% ggplot(aes(x=day)) +
#  geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
#  geom_line(aes(y=prediction, color=parent.pop)) +
#  facet_wrap(~parent.pop)

Min D1-13 (fit5)

all_wl2_size_for_models_temp_fit5 <- all_wl2_size_for_models_climate %>% drop_na(min_temp_d1_13)

temp_fit5 <- brm(formula=temp_f5, data=all_wl2_size_for_models_temp_fit5, prior=prior1_wl2, iter = 3800, cores=4)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG   -I"/usr/lib/R/site-library/Rcpp/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/unsupported"  -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/"  -I"/usr/lib/R/site-library/StanHeaders/include/"  -I"/usr/lib/R/site-library/RcppParallel/include/"  -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
##                  from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
##                  from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
##                  from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
##   679 | #include <cmath>
##       |          ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1
## Start sampling
summary(temp_fit5)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta) 
##          delta ~ 1
##          k ~ min_temp_d1_13
##          Hmin ~ (1 | parent.pop) + (1 | Genotype)
##          Hmax ~ (1 | parent.pop) + (1 | Genotype)
##    Data: all_wl2_size_for_models_temp_fit5 (Number of observations: 5403) 
##   Draws: 4 chains, each with iter = 3800; warmup = 1900; thin = 1;
##          total post-warmup draws = 7600
## 
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept)     1.77      0.05     1.68     1.87 1.00     3716     5460
## sd(Hmax_Intercept)     7.40      1.12     5.69     9.98 1.00     2465     3502
## 
## ~parent.pop (Number of levels: 22) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept)     1.51      0.26     1.08     2.08 1.00     3869     4435
## sd(Hmax_Intercept)    11.12      2.32     7.53    16.48 1.00     3028     3924
## 
## Regression Coefficients:
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## delta_Intercept      1.12      0.05     1.02     1.23 1.00     2682     4907
## k_Intercept          0.57      0.10     0.39     0.77 1.00     2386     3483
## k_min_temp_d1_13     0.01      0.00     0.00     0.01 1.00     6306     4789
## Hmin_Intercept       3.28      0.33     2.63     3.95 1.00     2372     3972
## Hmax_Intercept       8.92      2.57     4.25    14.59 1.00     1517     2538
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.76      0.01     0.74     0.77 1.00     5861     5796
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(temp_fit5, nvariables = 3, ask=FALSE)

#temp_fit5.predictions <- cbind(pred.df, prediction=predict(temp_fit5, newdata = pred.df)[,"Estimate"]) %>%
#  full_join(all_wl2_size_for_models, by=c("day", "parent.pop"))
#
#temp_fit5.predictions %>% ggplot(aes(x=day)) +
#  geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
#  geom_line(aes(y=prediction, color=parent.pop)) +
#  facet_wrap(~parent.pop)

Mean D1-13 (fit6)

all_wl2_size_for_models_temp_fit6 <- all_wl2_size_for_models_climate %>% drop_na(mean_temp_d1_13)

temp_fit6 <- brm(formula=temp_f6, data=all_wl2_size_for_models_temp_fit6, prior=prior1_wl2, iter = 3800, cores=4)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG   -I"/usr/lib/R/site-library/Rcpp/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/unsupported"  -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/"  -I"/usr/lib/R/site-library/StanHeaders/include/"  -I"/usr/lib/R/site-library/RcppParallel/include/"  -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
##                  from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
##                  from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
##                  from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
##   679 | #include <cmath>
##       |          ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1
## Start sampling
summary(temp_fit6)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta) 
##          delta ~ 1
##          k ~ mean_temp_d1_13
##          Hmin ~ (1 | parent.pop) + (1 | Genotype)
##          Hmax ~ (1 | parent.pop) + (1 | Genotype)
##    Data: all_wl2_size_for_models_temp_fit6 (Number of observations: 5403) 
##   Draws: 4 chains, each with iter = 3800; warmup = 1900; thin = 1;
##          total post-warmup draws = 7600
## 
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept)     1.78      0.05     1.68     1.87 1.00     4587     5397
## sd(Hmax_Intercept)     7.02      0.97     5.51     9.29 1.00     2045     3978
## 
## ~parent.pop (Number of levels: 22) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept)     1.49      0.26     1.08     2.07 1.00     3420     4655
## sd(Hmax_Intercept)    10.60      2.14     7.26    15.54 1.00     2800     4449
## 
## Regression Coefficients:
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## delta_Intercept       1.11      0.05     1.02     1.22 1.00     3031     4472
## k_Intercept           0.59      0.10     0.41     0.79 1.00     1907     3753
## k_mean_temp_d1_13     0.00      0.00     0.00     0.01 1.00     7555     4460
## Hmin_Intercept        3.27      0.33     2.64     3.94 1.00     2278     3673
## Hmax_Intercept        8.51      2.48     4.04    13.78 1.00     1544     2947
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.76      0.01     0.74     0.77 1.00     6782     5845
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(temp_fit6, nvariables = 3, ask=FALSE)

#temp_fit6.predictions <- cbind(pred.df, prediction=predict(temp_fit6, newdata = pred.df)[,"Estimate"]) %>%
#  full_join(all_wl2_size_for_models, by=c("day", "parent.pop"))
#
#temp_fit6.predictions %>% ggplot(aes(x=day)) +
#  geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
#  geom_line(aes(y=prediction, color=parent.pop)) +
#  facet_wrap(~parent.pop)

Max D1-13 (fit7)

all_wl2_size_for_models_temp_fit7 <- all_wl2_size_for_models_climate %>% drop_na(max_temp_d1_13)

temp_fit7 <- brm(formula=temp_f7, data=all_wl2_size_for_models_temp_fit7, prior=prior1_wl2, iter = 3800, cores=4)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG   -I"/usr/lib/R/site-library/Rcpp/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/unsupported"  -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/"  -I"/usr/lib/R/site-library/StanHeaders/include/"  -I"/usr/lib/R/site-library/RcppParallel/include/"  -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
##                  from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
##                  from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
##                  from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
##   679 | #include <cmath>
##       |          ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1
## Start sampling
summary(temp_fit7)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta) 
##          delta ~ 1
##          k ~ max_temp_d1_13
##          Hmin ~ (1 | parent.pop) + (1 | Genotype)
##          Hmax ~ (1 | parent.pop) + (1 | Genotype)
##    Data: all_wl2_size_for_models_temp_fit7 (Number of observations: 5403) 
##   Draws: 4 chains, each with iter = 3800; warmup = 1900; thin = 1;
##          total post-warmup draws = 7600
## 
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept)     1.79      0.05     1.69     1.88 1.00     4632     5990
## sd(Hmax_Intercept)     6.45      0.80     5.19     8.28 1.00     2471     3549
## 
## ~parent.pop (Number of levels: 22) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept)     1.44      0.26     1.03     2.03 1.00     3958     4880
## sd(Hmax_Intercept)     9.87      1.95     6.83    14.47 1.00     3535     5059
## 
## Regression Coefficients:
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## delta_Intercept      1.06      0.04     0.98     1.14 1.00     3159     4830
## k_Intercept          0.66      0.10     0.47     0.86 1.00     2341     3512
## k_max_temp_d1_13     0.00      0.00     0.00     0.00 1.00     5504     4237
## Hmin_Intercept       3.25      0.32     2.64     3.89 1.00     2790     3996
## Hmax_Intercept       7.91      2.30     3.54    12.76 1.00     1490     2808
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.76      0.01     0.74     0.77 1.00     6152     6133
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(temp_fit7, nvariables = 3, ask=FALSE)

#temp_fit7.predictions <- cbind(pred.df, prediction=predict(temp_fit7, newdata = pred.df)[,"Estimate"]) %>%
#  full_join(all_wl2_size_for_models, by=c("day", "parent.pop"))
#
#temp_fit7.predictions %>% ggplot(aes(x=day)) +
#  geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
#  geom_line(aes(y=prediction, color=parent.pop)) +
#  facet_wrap(~parent.pop)
save.image("../output/BQC_WL2_Growth_Climate.Rdata")

Moisture Models

moisture_f2 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
                  delta ~ 1,
                  k ~ s_moisture_1_7, 
                  Hmin + Hmax ~ (1|parent.pop) + (1|Genotype),
                  nl=TRUE)

moisture_f3 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
                  delta ~ 1,
                  k ~ s_moisture_1_13,
                  Hmin + Hmax ~ (1|parent.pop) ,
                  nl=TRUE)

Moisture Model Fits

Min D1-7 (fit2)

all_wl2_size_for_models_moisture_fit2 <- all_wl2_size_for_models_climate %>% drop_na(s_moisture_1_7)

moisture_fit2 <- brm(formula=moisture_f2, data=all_wl2_size_for_models_moisture_fit2, prior=prior1_wl2, iter = 3800, cores=4)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG   -I"/usr/lib/R/site-library/Rcpp/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/unsupported"  -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/"  -I"/usr/lib/R/site-library/StanHeaders/include/"  -I"/usr/lib/R/site-library/RcppParallel/include/"  -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
##                  from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
##                  from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
##                  from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
##   679 | #include <cmath>
##       |          ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1
## Start sampling
summary(moisture_fit2)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta) 
##          delta ~ 1
##          k ~ s_moisture_1_7
##          Hmin ~ (1 | parent.pop) + (1 | Genotype)
##          Hmax ~ (1 | parent.pop) + (1 | Genotype)
##    Data: all_wl2_size_for_models_moisture_fit2 (Number of observations: 6038) 
##   Draws: 4 chains, each with iter = 3800; warmup = 1900; thin = 1;
##          total post-warmup draws = 7600
## 
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept)     1.86      0.05     1.76     1.96 1.00     4761     4988
## sd(Hmax_Intercept)     5.20      0.41     4.51     6.14 1.00     2774     3997
## 
## ~parent.pop (Number of levels: 22) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept)     1.39      0.25     0.99     1.94 1.00     4104     4462
## sd(Hmax_Intercept)     8.44      1.43     6.08    11.71 1.00     4223     5303
## 
## Regression Coefficients:
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## delta_Intercept      1.03      0.04     0.95     1.12 1.00     2257     4565
## k_Intercept          0.81      0.09     0.62     0.99 1.00     2609     3720
## k_s_moisture_1_7     0.13      0.08     0.02     0.31 1.00     8524     4542
## Hmin_Intercept       3.25      0.31     2.64     3.89 1.00     2939     4300
## Hmax_Intercept       7.08      1.89     3.48    10.95 1.00     1889     3346
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.84      0.01     0.82     0.86 1.00     5407     5643
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(moisture_fit2, nvariables = 3, ask=FALSE)

#moisture_fit2.predictions <- cbind(pred.df, prediction=predict(moisture_fit2, newdata = pred.df)[,"Estimate"]) %>%
#  full_join(all_wl2_size_for_models, by=c("day", "parent.pop"))
#
#moisture_fit2.predictions %>% ggplot(aes(x=day)) +
#  geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
#  geom_line(aes(y=prediction, color=parent.pop)) +
#  facet_wrap(~parent.pop)

Mean D1-13 (fit3)

all_wl2_size_for_models_moisture_fit3 <- all_wl2_size_for_models_climate %>% drop_na(s_moisture_1_13)

moisture_fit3 <- brm(formula=moisture_f3, data=all_wl2_size_for_models_moisture_fit3, prior=prior1_wl2, iter = 3800, cores=4)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG   -I"/usr/lib/R/site-library/Rcpp/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/"  -I"/usr/lib/R/site-library/RcppEigen/include/unsupported"  -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/"  -I"/usr/lib/R/site-library/StanHeaders/include/"  -I"/usr/lib/R/site-library/RcppParallel/include/"  -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
##                  from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
##                  from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
##                  from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
##   679 | #include <cmath>
##       |          ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1
## Start sampling
summary(moisture_fit3)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta) 
##          delta ~ 1
##          k ~ s_moisture_1_13
##          Hmin ~ (1 | parent.pop)
##          Hmax ~ (1 | parent.pop)
##    Data: all_wl2_size_for_models_moisture_fit3 (Number of observations: 6038) 
##   Draws: 4 chains, each with iter = 3800; warmup = 1900; thin = 1;
##          total post-warmup draws = 7600
## 
## Multilevel Hyperparameters:
## ~parent.pop (Number of levels: 22) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept)     0.97      0.24     0.59     1.52 1.00     1946     2974
## sd(Hmax_Intercept)     5.79      1.01     4.22     8.06 1.00     2056     2872
## 
## Regression Coefficients:
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## delta_Intercept       0.95      0.13     0.74     1.23 1.00     2043     2328
## k_Intercept           1.93      0.28     1.31     2.41 1.00     5820     3273
## k_s_moisture_1_13     0.70      0.45     0.10     1.80 1.00     9792     4541
## Hmin_Intercept        2.89      0.28     2.37     3.47 1.00     2176     2590
## Hmax_Intercept        5.72      1.26     3.21     8.35 1.01      732     1282
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     2.49      0.02     2.44     2.53 1.00    13345     5406
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(moisture_fit3, nvariables = 3, ask=FALSE)

#moisture_fit3.predictions <- cbind(pred.df, prediction=predict(moisture_fit3, newdata = pred.df)[,"Estimate"]) %>%
#  full_join(all_wl2_size_for_models, by=c("day", "parent.pop"))#

#moisture_fit3.predictions %>% ggplot(aes(x=day)) +
#  geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
#  geom_line(aes(y=prediction, color=parent.pop)) +
#  facet_wrap(~parent.pop)

Temp Model Comparison

#base_fit1 <- add_criterion(base_fit1, "loo")
#base_fit2 <- add_criterion(base_fit2, "loo")
temp_fit2 <- add_criterion(temp_fit2, "loo")
## Warning: Found 331 observations with a pareto_k > 0.7 in model 'temp_fit2'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
temp_fit3 <- add_criterion(temp_fit3, "loo")
## Warning: Found 319 observations with a pareto_k > 0.7 in model 'temp_fit3'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
temp_fit4 <- add_criterion(temp_fit4, "loo")
## Warning: Found 324 observations with a pareto_k > 0.7 in model 'temp_fit4'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
temp_fit5 <- add_criterion(temp_fit5, "loo") 
## Warning: Found 328 observations with a pareto_k > 0.7 in model 'temp_fit5'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
temp_fit6 <- add_criterion(temp_fit6, "loo") 
## Warning: Found 340 observations with a pareto_k > 0.7 in model 'temp_fit6'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
temp_fit7 <- add_criterion(temp_fit7, "loo") 
## Warning: Found 324 observations with a pareto_k > 0.7 in model 'temp_fit7'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
loo_compare(temp_fit2, temp_fit3, temp_fit4, temp_fit5, temp_fit6, temp_fit7)
##           elpd_diff se_diff
## temp_fit4   0.0       0.0  
## temp_fit3  -1.7       7.2  
## temp_fit7  -5.2      10.1  
## temp_fit6  -5.6      10.3  
## temp_fit5  -6.7      13.3  
## temp_fit2 -13.0      12.3

Moisture Model Comparison

#base_fit1 <- add_criterion(base_fit1, "loo")
#base_fit2 <- add_criterion(base_fit2, "loo")
moisture_fit2 <- add_criterion(moisture_fit2, "loo")
## Warning: Found 291 observations with a pareto_k > 0.7 in model 'moisture_fit2'.
## We recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
moisture_fit3 <- add_criterion(moisture_fit3, "loo")
loo_compare(moisture_fit2, moisture_fit3)
##               elpd_diff se_diff
## moisture_fit2     0.0       0.0
## moisture_fit3 -5670.1     179.1
save.image("../output/BQC_WL2_Growth_Climate.Rdata")